A polarizable interatomic force field for Ti02 parameterized using density 
functional theory 
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We report a classical interatomic force field for Ti02, which has been parameterized using density 
functional theory forces, energies, and stresses in the rutile crystal structure. The reliability of 
this new classical potential is tested by evaluating the structural properties, equation of state, 
phonon properties, thermal expansion, and some thermodynamic quantities such as entropy, free 
energy, and specific heat under constant volume. The good agreement of our results with ab 
initio calculations and with experimental data, indicates that our force-field describes the atomic 
interactions of Ti02 in the rutile structure very well. The force field can also describe the structures 
of the brookite and anatase crystals with good accuracy. 



PACS numbers: 34.20. Cf, 31.15.es, 63.20.-e, 64.70. kp 



I. INTRODUCTION 

Recently, a bipolar switching phenomenon in TiC>2 has 
been observed yj-y) and much work has been done to de- 
sign superior TiC>2-based resistive random access mem- 
ory (RRAM) devices. This, together with many other 
important applications of TiC>2, such as white pigment 
for paints, sunscreen, high-efficiency solar cells (|J), and 
photosplitting of water to hydrogen has been stim- 

ulating a great deal of research interest in atomistic sim- 
ulations of TiC>2. Accurate atomistic simulations would 
provide detailed information at the atomic level that 
might answer some important open questions, such as 
the mechanism of bipolar switching. 

Force fields are the key to accurate atomistic simula- 
tions. In the last twen ty y ears, many force fields have 
been built for TiC>2 (|7H14). Among these, the partial- 
charge rigid-ion model of Matsui and Akaogi (MA model) 
(0), and the variable-charge Morse stretch (MS-Q) inter- 
atomic potential of Swamy|8|), have been the most suc- 
cessful and popular. Both the MA and MS-Q models can 
describe a series of crystal structures quite well. A com- 
mon feature of these force fields is that they were adapted 
to reproduce experimental bulk properties, such as lat- 
tice constants, cohesive energies, bulk moduli and elastic 
constants. An alternative parameterization method con- 
sists of fitting force fields to forces, energy differences, 
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and stresses extracted from ab initio calculations (fl5l- 
0j|). To our knowledge, no such ab initio force field for 
TiC>2 yet exists. 

Having a potential that can accurately model the vi- 
brational properties of TiC>2 as well as its structures is 
very important as they determine its finite-temperature 
macroscopic properties. The phonon dispersion curves 
and density of states for rutile TiC>2 have been measured 
by the coherent inelastic scattering of thermal neutrons 
along principal symmetry directions of the Brillouin zone 
(fl9l). The T-point phonons have been measured by Ra- 
man (pol ) and infrared spectroscopy. In addition to 
these experimental determinations,phonon frequencies 
have been calculated ab initio \XM2M ) . However, little 
attention has been paid to the determination of phonon 
properties from a force field for TiC>2, even though it is 
a very demanding test of a force field's accuracy. As has 
been demonstrated for MgQ (fl7l: [25h . ab initio-b&sed clas- 
sical force fields have the potential to give a very good 
description of phonon properties. TiC>2 has several com- 
peting lattice structures and is a more complex oxide 
than MgO. It is rather interesting, therefore, to see how 
well such a force field can predict its vibrational proper- 
ties. 

In this work, we present a polarizable classical inter- 
atomic force field for Ti02, parameterized from ab ini- 
tio calculations. We test its accuracy by comparing its 
predictions of structural properties, equations of state, 
phonon frequencies, and the temperature dependences of 
lattice parameters, entropy, free energy, and specific heat 
under constant volume in the rutile phase with experi- 
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mental data and with ab initio calculations. We also cal- 
culate the structures and relative energies of the brookite 
and anatase phases to get an indication of the force field's 
transferability. 



II. FORCE FIELD 

The potential energy function U = U SI + U cs that 
we use to describe interactions between ions consists 
of a pairwise term U sr describing short-range non- 
electrostatic interactions and an electrostatic term U cs 
describing dipole induction and the electrostatic interac- 
tions of the charges and induced dipole moments of the 
ions. 

For the short-range interaction energy we use the pair- 
wise Morse-stretch form 
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where ry = |r,- — rj\ is the distance between nearby ions 
i and j and , 7^ , and are parameters that are spe- 
cific to the pair of species of ions i and j. This potential 
is truncated at a radius of 18.0 a.u. 

The total electrostatic contribution to the energy of the 
system includes charge-charge, charge-dipole, and dipolc- 
dipole interaction terms: 
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where qi is the charge of ion i, pj is the 7 th cartesian 
component of the dipole moment of ion «, a, is its polar- 
isability, and V 7 = d/drj^. The final term on the right 
hand side of this equation is a sum of the self energies of 
the ions. The self energy of an ion is the internal energy 
cost of inducing a dipole on it. 

There are two mechanisms for the induction of ionic 
dipole moments. The first is via the short-range repulsive 
forces between ions which can distort an ion's electron 
cloud thereby giving it a dipole moment. The second is 
the induction of a dipole moment on an ion by the electric 
field arising from the charges and the dipole moments of 
all other ions. 

To model the induction by short-range re puls ive forces 
we follow the approach of Madden et aZ. <|26l; l27h ■ In their 
approach, the contribution of the short-range forces to 
dipole moments is given by 
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and bij and are model parameters. 

The electrostatically-induced dipole moments are more 
difficult to calculate because they are all interdependent 
due to their contributions to, and dependences on, the 
electric field. At every time step we find the dipole 
moments of the ions by iterating to self-consistency the 
equation 
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where p™ is the dipole moment on ion i at iteration m 
of the self-consistent cycle; E(r,;) is the electric field at 
position rj. 

We use Ewald summation to calculate electrostatic en- 
ergies, forces, stresses, and electric fields. 

Our model depends on a set of parameters {r] n } com- 
prising the Morse-Stretch parameters D, 7, and r , the 
parameters b and c of the short-range dipole model, and 
the charges q Ti , q a and polarisabilities a Ti , a Q of the ions. 
In section HITl we describe how these parameters are de- 
termined. 



III. POTENTIAL PARAMETERIZATION 

A. Ab initio molecular dynamics 

Total energies, stresses, and forces are computed from 
ab initio simulations based on density functional the- 
ory (DFT). Ab initio molecular dynamics (MD) simu- 
lations employing the projector augmented wave method 
(PAW) were performed in the NVT ensemble under pe- 
riodic boundary conditions using the VASP simulation 
package (JH^. Only the Ti (3d,4s) and O (2s, 2p) elec- 
tron states were treated as valence states, however, tests 
in which the Ti (3s, 3p) semicore states were included 
yielded similar results for the energy differences between 
crystalline phases and the phonon frequencies of rutile 
Ti02- The local-density approximation (LDA) to the 
exchange-correlation energy was used. The velocity Ver- 
let algorithm (|3ll ) with a time step of 1 fs was adopted 
to solve the equations of motion, and temperature was 
controlled using a Nose-Hoover thermostat (|32l). We used 
a 2 x 2 x 3 superccll that contained 72 atoms and only the 
T-point was used to sample the Brillouin zone. A kinetic 
energy cutoff of 180 eV was used in the plane-wave ex- 
pansion of the wave functions during the MD simulations. 
The goal of these ab initio MD simulations was simply to 
produce some reasonable atomic configurations to serve 
as representative snapshots at finite temperature. These 
configurations were then used to perform higher-precision 
DFT calculations of forces, stresses, and energies for use 
in the force-field parameterization. 



B. Fitting procedure for the potential 

We fit our potential to DFT calculations in the ru- 
tile crystal structure, however, to broaden the range of 
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different environments for the ions, two temperatures are 
considered in the potential fitting: 300K and 1000K. Fur- 
thermore, at each temperature, configurations at more 
than one volume were chosen. For 300K eight configu- 
rations were used. Six of these configurations are taken 
from an ab initio MD simulation at the zero temperature 
equilibrium volume. The other two were obtained from 
simulations in which the volume was expanded by 3%. At 
1000K, we used 3 sets of data in total, each containing 
four configurations. These three sets are obtained from 
MD simulations at the zero-temperature equilibrium vol- 
ume, and at volumes expanded by 3% and by 9%. To 
minimize correlations between configurations, successive 
snapshots were separated from one another by 2 ps of ab 
initio MD. For each snapshot, we used VASP to calculate 
the forces, energies and stresses. To obtain these quanti- 
ties with high precision, the Brillouin zone was sampled 
using a 2 x 2 x 2 Monkhorst-Pack k— point mesh and, to 
converge the forces and the stress tensors properly, a very 
high wavefunction cutoff energy of 1500 eV was used. 

The potential parameters were obtained by minimizing 
the cost function 

r({?7 n }) = w f AF + w s AS + w e AE (6) 

with respect to the parameters {r] n } where 
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F^ i is the a-component of the force on ion i calculated 
with the classical potential; F^ i is the a-component of 

the force on ion i calculated ab initio; S^f is the stress 
tensor component calculated with the classical potential; 
S"f is the stress tensor component calculated ab initio; 
N c denotes the number of snapshot configurations used in 
the fit; B is the bulk modulus; Uj^ is the potential energy 
of configuration k for the classical potential, and is its 
energy calculated ab initio. Wf, w s , and w e are weights 
that quantify the relative importances to the fitting cost 
function of forces, stresses, and energies, respectively. In 
this work, Wf = 1.0, w s = 0.5, and w e = 0.01. These 
numbers are somewhat arbitrary as long as the energy 
differences are given a relatively small weight because we 
only have one energy per configuration. We tried differ- 
ent weights and found no significant difference in the final 
values of AF, AS, and AE. Minimization of r({77„}) 
is performed by a combination of simulated annealing 
((331 ) and Powell minimization (|34h . The former is used 
to provide an initial parameter set which brings the cost 
function to a basin in the surface defined by r({ri n }) in 



77— space. Minimization is then completed using Powell's 
method. 

We assume that DFT provides accurate atomic forces 
and energies. Therefore if, at the global minimum of T 
in 77— space, AF remains large it means that the poten- 
tial model does not contain the ingredients necessary to 
fit the DFT potential energy surface closely. The model 
is unphysical or incomplete and the resulting force-field 
is unlikely to produce results that are in good agreement 
with experiment. Agreement with experiment can still be 
found by fitting to experimental data, but when agree- 
ment is not attributable to accurate interatomic forces, 
it is difficult to place confidence in the force-field for cal- 
culations that can not be directly checked by comparison 
with experiment. However, when reliable experimental 
data is available, simulations are of limited value anyway. 
For this reason we deem it important for a force-field to 
be able to fit DFT forces closely and the closeness of the 
fit achieved in the parameterisation is itself a test of the 
quality of the potential. 



IV. RESULTS 

Following the fitting procedure described in the previ- 
ous section, we obtained the force-field parameters for ru- 
tile Ti0 2 , as listed in table HI In this fitting, AF = 0.202, 
AS = 0.009, and AF = 0.046. Put another way, the 
root-mean-squared error in the forces is 20.2% of the 
root-mean-squared force. 

As a comparison we have checked how close the MA 
model's forces are to the ab initio ones. We find that 
the root-mean-squared difference between the MA forces 
and the DFT forces is ~ 88% of the root-mean-squared 
DFT forces (AF w 0.88, AS 0.02, AF « 0.63). The 
MA model uses the pairwise Born-Mayer potential form 
and by minimising T with respect to the parameters of 
the Born-Mayer potential we have been able to achieve 
a best fit of AF « 0.48, AS « 0.01, AF w 0.10. How- 
ever, this resulted in a disimprovement in the description 
of the three crystal structures with respect to the MA 
potential. As one might expect, and as was found previ- 
ously for silica lflTjh . when the fit to ab initio calculations 
is poor, the ability of a potential to reproduce experi- 
mental results does not correlate with the quality of the 
fit. Our tests strongly suggest that the Born-Mayer form 
provides a poor microscopic description of atomic forces 
in TiC>2. 



A. Crystal structures 

In order to check the reliability and transferability of 
the new force field, we computed the structural properties 
of three different crystalline polymorphs of Ti02, namely 
the rutile, anatase, and brookite structures. Structural 
optimizations for these three crystals were performed via 
the steepest-descent method at zero pressure. The re- 
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Parameters 


O-O 


Ti-0 


Ti-Ti 


D 


1.5865 x 10" 3 


1.7668 x 10 -3 


5.0276 x 10" 3 


7 


9.6370 


12.2332 


5.9069 


r° 


6.5150 


4.7678 


8.1380 


b 


1.5368 


4.8187 


1.8246 


c 


2.9216 


-122.1489 


-0.7918 


a 


2.9314 




10.2739 


q 


-1.1045 




2.209 



TABLE I The fitted interatomic force-field parameters (in 
atomic units). 
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FIG. 1 Volume dependence of the energy of TiC>2 in the ru- 
tile, brookite, and anatase crystal structures for our potential 
(red lines) and from LDA ab initio calculations. Energies are 
expressed relative to the equilibrium energy of the most stable 
phase. 



suiting lattice parameters are presented in Table |TT] and 
compared to the ab initio data calculated with the VASP 
package and to experimental data (|35l437t ) . both at room 
temperature and (for rutile and anatase) at 15 K. Both 
ab initio calculations and the new force field reproduce 
experimental data quite well and the accuracy of the 
force field is comparable to that of DFT. We stress that 
the force field in this work is only fit to ab initio data 
from calculations in the rutile crystal structure. There- 
fore, its ability to reproduce the structures of anatase 
and brookite is an indication that it may be reasonably 
transferable between different structures. However, apart 
from the handful of structural and energetic properties 
presented in this paper, we have not tested the force field 
in these phases. 



B. Equations of state 

The energy as a function of volume for all three phases 
has been calculated with our potential and the results 
are illustrated in Fig. [TJ Structural relaxations were 
performed at each volume. The results in red are from 



the new force field, while those in blue are from our 
density functional theory calculations within the local 
density approximation (DFT-LDA). In each case, we 
set the energy of the most stable phase in equilibrium 
to zero. We can see that the DFT-LDA calculations 
and the new force field give similar equilibrium vol- 
umes for the three crystals (Table HTH) . However, it is 
clear from Fig. Q] and Table IIV1 that the energy se- 
quences obtained from DFT-LDA calculations and the 
new force field are different and, furthermore, that the 
the magnitudes of the energy differences between the 
three structures differ substantially between DFT-LDA 
and our force field. In DFT-LDA calculations, the en- 
ergy differences between the three structures are small, 
and -EWookitc < ^anatase < futile- However, for the new 
potential, the rutile structure is the ground state, and 
-Erutiic < ^brookite < ^anatase- While this energy order- 
ing is consistent with the experimental results (|38l ; 1391 ) 
the magnitudes of the energy differences are at least an 
order of magnitude too large. The DFT-LDA estimates 
of the magnitudes of the energy differences, both here and 
in previous work(|40|). are in much better agreement with 
experiment, despite the fact that they predict the wrong 
sequence of structures. We have also calculated the equi- 
librium energy differences within a generalized gradient 
approximation (GGA) (j4l[ ) to the exchange-correlation 
energy. The GGA differences in energy between anatase 
and rutile and between brookite and rutile are also pre- 
sented in Table IIV1 The magnitudes of these energy dif- 
ferences are larger than those of LDA, but still a factor 
of ~ 5 smaller than we find with our force field. We 
are unaware of compelling reasons why LDA should out- 
perform GGA for Ti02- Therefore, it seems reasonable 
to consider differences between LDA and GGA as lower 
bounds on the error bars associated with our approxima- 
tion to the exchange-correlation energy. The free-energy 
differences calculated with the MA model at room tem- 
perature are similar in magnitude (but slightly smaller) 
to the K energy differences calculated with our force 
field and the MA model's energy ordering of the three 
crystal structures is also consistent with experiment. Al- 
though our force field produces energy differences that 
are much too large, it is at least gratifying that the crys- 
tal structures have the correct energy sequence. 

The relationship between energy and volume can be 
described by the Birch-Murnaghan equation of state (fi3 ) 



E(V)-E = 
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where V is the volume, Vq the equilibrium volume, Eg 
the equilibrium energy, B the bulk modulus, and Bq the 
pressure derivative of the bulk modulus. The parameters 
for the equation of state are listed in Table IIII| together 
with experimental data (@; H^; [44]) . The results from the 
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~0K 




9.1197 


9.095 




5.415 


5.399 




5.103 


5.131 


~ 300 K 


9.174(37) 






5.449(37) 






5.138(37) 







TABLE II Lattice parameters for three optimised structures of Ti02, together with experimental data and results from ab initio 
calculations. DFT-LDA represents ab initio calculations with the VASP package within the local density approximation, and 
POT the new force field, obtained by fitting to ab initio data in the rutile structure. Data at two temperatures are presented, 
where available: at close to K (15 K in the experiment of Ref. [36j, K for simulation) and at room temperature. 





Vo/TiOa 
DFT-LDA 


(A 3 ) 

POT 


Exp. 


Bo (GPa) 
DFT-LDA 


POT 


DFT-LDA 


POT 


Rutile 


30.62 


30.46 


211± 7(43) 


252 


241 


5.37 


5.07 








230± 20(44) 










Anatase 


33.58 


33.37 


178(8) 


188 


171 


1.75 


2.34 


Brookite 


31.55 


31.50 


no data 


230 


196 


3.87 


3.14 



TABLE III Birch-Murnaghan EOS parameters for different structures of Ti02, together with the experimental data available. 
Experimental data refers to room temperature, while theoretical values are obtained at K. 



ab initio calculations, the new classical potential, and the 
experimental data are all in good agreement, however, it 
is important to note that the simulation results refer to 
zero temperature while the experiments were performed 
at room temperature. 



C. Vibrational properties of rutile 

To have a further test of the new force field, the 
phonons of rutile Ti02 are calculated with the small 
displacement method implemented by Togo's FROPHO 
package (fill). Figure [2^ a) compares phonons calculated 
with our classical potential (black dots) with those cal- 
culated from DFT-LDA (green triangles) and, for some 
branches and symmetry directions, with inelastic neutron 
data at room temperature taken from reference 19. Fig- 
ure J^b) shows the same ab initio and experimental data, 
but the black dots are now from calculations using the 
MA model (|7|). We include in both these plots only those 
frequencies calculated for phonons whose wavelengths are 
commensurate with our simulation cells, i.e., no interpo- 
lation has been performed to infer frequencies at differ- 
ent wavevectors. Simulations with the classical poten- 
tials were performed with three different simulation su- 
percells whose sizes, in multiples of the primitive 6-atom 
unit cell, were 14 x 14 x 14 (16464 atoms), 12 x 12 x 12 



(10368 atoms), and 10 x 10 x 10 (6000 atoms). Our DFT 
calculations were performed in a 4 x 4 x 4 (384 atoms) 
supercell. 



In our calculations we have corrected longitudinal op- 
tical frequencies in the long wavelength limit (i.e. k — » 0) 
to account for the effects of the long-range electric fields 
that these modes induce (fdH lifjh. 



The first thing to note from Figures [Ha) and (Hb) is 
that both classical potentials struggle to reproduce the 
dispersion of the acoustic modes close to the zone bound- 
aries. This can be seen most clearly along the T — > A and 
r — > Z symmetry directions. However, what is also clear 
is that our force field greatly improves on the MA model 
in this respect. For the high-frequency modes, between 
~ 15 THz and ~ 25 THz, our model is in very good agree- 
ment with experiment and with DFT. The MA model's 
description of the phonon spectrum at these frequencies 
is very poor by comparison. Overall, it is clear that our 
potential provides a description of the phonon spectrum 
that is much better than that provided by the MA model. 
It is also clear, however, that there is room for improve- 
ment, presumably by improving the functional form of 
the potential. 
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DFT-LDA 


DFT-GGA 


POT 


MA 


E ana tasc E ru tjl e 


35.0(38), 53.8(39) 


-12.0 


-81.2 


425.9 


301.6 


Ebrookite~E ru tile 


7.7(38) 


-17.5 


-39.9 


212.6 


179.3 



TABLE IV Energy differences in meV/Ti02 between the anatase and brookite crystal structures and the rutile crystal structure. 
The experiments of reference [3^ and the calculations with the MA model(0) were performed at room temperature, while the 
more recent experiments of reference [33 were performed at 971 K. The experiments measured enthalpies of transformation 
whereas all the DFT and force-field calculations that we report are energy differences at K. 
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FIG. 2 Phonon dispersion curves of Ti02 in the rutile crystal 
structure from different methods. (a)Results from our new 
force field (black dots) compared to experiment lfl9h (red cir- 
cles) and ab initio (DFT-LDA) calculations (green triangles). 
(b)Results from the MA force field (black dots) compared to 
experiment (jl9|) (red circles) and ab initio (DFT-LDA) calcu- 
lations (green triangles). 



D. "Electronic" properties 

To correct the LO frequencies at the T-point, it has 
been necessary to calculate the high-frequency dielec- 
tric tensor Eqq, which is diagonal, and the Born effective 
charge tensor Z*. We have calculated £oo and Z* directly 
by calculating the polarization responses to small elec- 
tric fields (|4!^) and to small displacements of the atoms, 
respectively. In the rutile crystal structure we find the 
components of perpendicular and parallel to the c- 
axis to be e± — 2.58 and £|| = 2.75, respectively, which 
is in poor agreement with those calculated ab initiol^di) 
(e± = 7.49, £|| = 8.57). The Born effective charge ten- 
sor of rutile has only three independent non-zero compo- 
nents for each atom. We calculate them to be Z* x (Ti) = 
Z*,.(Ti) = 3.428, Z* y (Ti) = Z* (Ti) = 0.422, Z* zz (Ti) 



yy 

3.789 and 

Z* VX (0) = 



z* x (of y = z; y (d) = -1.714, V*(o) = 

-1.044, Z* z (0) = -1.895. Agaimthis is 
in poor agreement with the ab initio results (50,) which 
are Z* x (Ti) = Z* yy (Ti) = 6.36, Z* xy (Ti) = Z* yx (Ti) = 
1.00, Z* zz (Ti) = 7.52 and Z* x (0) = Z* yy (0) = -3.18, 
Z* y (0) = Z; x (0) = -1.81, Z* z (0) = -3.76. Our po- 
tential is not intended to give a good description of the 
electronic properties of the system, so this disagreement 
is no cause for concern. We intend, instead, an accu- 
rate representation of the interatomic forces. Both Z* 
and £oo are derivatives of the polarization field, and are 
therefore only meaningful on length scales that are large 
compared to the primitive unit cell of the crystal. At 
large distances, forces between atoms (labeled 1 and 2) 
are proportional to Z^Z^/Soo, and so this is the relevant 
quantity if one is interested assessing the quality of these 
response functions' contributions to interatomic forces. 
Therefore, we look at the components of the screened 



effective charge tensor 



The symmetry of 



the crystal is such that, for each atomic species, there 
are only three non-unique non-zero components of the 
screened effective charge tensor: Z xx /^/e± = Z* y /^/e±, 
Z xy/Ve± = Z yx/Ve±> and Z lz/VH- 

We have calculated these quantities, which determine 
long-range forces, with our force field and, in Table [V] 
we compare our results to those obtained from DFT. 
The agreement between our potential and DFT is good, 
therefore our potential should describe long-range forces 
reasonably well. 

The best-fit ionic polarisabilities ( ao = 2.9 a.u., 
aTi = 10.3 a.u.) (Table HJ) provided by our parame- 
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Ti 


DFT 





Ti 


POT 


O 




/eX = 




2.32 




-1.16 


2.13 




-1.07 




/ex = 


Zy X / \fS~L 


0.37 




-0.66 


0.26 




-0.65 






m 


2.57 




-1.28 


2.28 




-1.14 



TABLE V Components of the screened effective charge tensor calculated with density functional theory (DFT) and with our 
new potential (POT) 



terization process are unusual and worthy of comment. 
The polarisability of the oxygen ion is much lower than 
has been found by force fitting for other oxides such 
as MgO(|l7|) and SiQ^(|16l). or from electronic structure 
simulations(25). Typically ao is between 5 and 15 a.u.. 
However, most surprising is the large polarizability of the 
Ti ions. One would expect a small polarizability for a 
cation, particularly one with such a high positive charge. 

By making atoms polarisable, we include the response 
of electrons in a phenomenological way. However, we 
only include one response mechanism. Others, such as 
higher-order polarisabilities and charge-transfer between 
ions, are certain to exist to some degree. Furthermore, we 
are using a fully-ionic model while covalent effects might 
be important. It seems likely that including one or more 
of these effects would allow our parameterizer to achieve 
a closer fit to the ab initio forces. We speculate that the 
strange polarisabilities are an effort by our parameteri- 
zation program to compensate for those other electronic 
effects which are not present in the mathematical form 
of the potential. It is remarkable, and further evidence 
of the power of the force-fitting approach, that the pa- 
ramcteriser succeeds in finding values for the charges and 
polarisabilities of the atoms that perform so well. 



E. Thermal expansion 

Figure 3 compares the dependence of the lattice pa- 
rameters of the rutile structure on temperature with 
experiment (|37l ; 48). Once again, the agreement is very 
good, indicating that the anharmonicity of the potential 
energy surface is well reproduced by our force field. The 
results reported in figure 3 were from long (> 100 ps) MD 
simulations of a 6 x 6 x 8 supercell, which contains 1728 
atoms. A time step of 0.723 fs was used and tempera- 
ture was controlled with a Nose- Hoover thermostat ((H). 
While atoms moved according to Newton's equations 
(using a Verlet algorithm), steepest descent was per- 
formed on the cell degrees of freedom. We verified that 
this approach gave almost identical results to perform- 
ing Parrinello- Rahman constant pressure simulations (f5lh 
and taking averages over the trajectory of the lattice con- 
stants. However, the latter simulations are more time 
consuming and more difficult to control. We caution that 
our simulations treat the ions purely classically and that 
quantum effects would flatten out these curves at low 
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Temperature [Kelvin] 
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Temperature [Kelvin] 



FIG. 3 Temperature dependence of the lengths of the ru- 
tile lattice vectors, expressed as percentage differences with 
respect to their room temperature values, i.e. Aa(T) — 
a(T) — a(300 K). Boxes and filled triangles refer to the ex- 
perimental data of Refs. 48 and respectively. Blue dots 
are the results of MD simulations with our new force field. 



temperatures. 



F. Thermodynamic properties 

Thermodynamic properties of rutile Ti02 , such as free 
energy, entropy, and phonon specific heat under constant 
volume were evaluated in the quasi-harmonic approxima- 
tion from the phonons using the FROPHO package (|45l ) 
and the results are given in Figure 2] We compare the 
results from DFT calculations of the phonons and calcu- 
lations of the phonons with our force field. In both cases, 
phonons were calculated in 4 x 4 x 4 supercell. As illus- 
trated in Figure HJa), the temperature dependence of the 
free energy calculated from the new force field is almost 
the same as that from our ab initio calculation. An en- 
largement of the free energy curve is also superimposed 
in Figure |U a) in order to see the difference between these 
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two results more clearly. The entropy of rutile Ti02 is 
given in Figure |4]Jb). Again, the new force field repro- 
duces the ab initio calculation extremely well. 

FigureHJc) shows that the specific heat under constant 
volume (C v ) from the new classical force field is also in ex- 
cellent agreement with our ab initio calculations. Within 
the Debye model C V (T) can be expressed as 



r, = <>A7,, ; ( I- 



t d /t 



(e* - 1) 



(11) 



where N is the number of atoms per cell (6 for rutile 
Ti02), fcs is the Boltzmann constant, and To is the De- 
bye temperature. Using this expression, the specific heat 
at To can be computed as 



C v (T = T D )=9Nk B 



1 



T dx = 17.1312fc B 



(12) 

By comparing this value with our calculation, we can 
obtain the Debye temperature To. We find that our ab 
initio calculations predict To = 781.6 K while our new 
force field predicts Td = 786.1 K. These values are very 
close to the experimental Debye temperature of 778.3 K 
d2£ 



V. CONCLUSIONS 

A new classical force field for TiC>2 has been pro- 
posed that has been parameterized using forces, stresses, 
and energies extracted from ab initio molecular dynam- 
ics simulations. The potential can not only predict the 
structural properties of Ti0 2 in three crystal structures, 
namely, rutile, brookite, and anatase, but also gives the 
right ground state and the correct energy sequence of 
these three structures and predicts their bulk moduli with 
an accuracy comparable to DFT. While the force field has 
been constructed with rutile in mind, and fit to ab initio 
data of the rutile structure, it seems likely to be reason- 
ably transferable to different structures. However, the 
most obvious deficiency of the force-field is that it over- 
estimates energy differences between these polymorphs. 
This, and the unphysical polarisabilities of the ions, sug- 
gest that our functional form may lack ingredients nec- 
essary to describe some electronic effects that can play 
important roles in the energetics of TiC>2. Quadrupole 
polarisation of ions and dynamical charge transfer be- 
tween ions are two effects that may have an important 
role to play. 

Concerning the vibrational properties, the new poten- 
tial can reproduce most features of ab initio and exper- 
imental results, and is a substantial improvement over 
the widely- used MA model for TiC>2. Furthermore, the 
potential provides an excellent description of the tem- 
perature dependence of the rutile lattice constants. The 
good description of vibrational properties from the pro- 
posed force field is further confirmed by thermodynamic 
quantities calculated from the phonon properties. The 
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FIG. 4 Thermodynamic properties of rutile TiC>2 derived 
from phonon properties: (a) free energy F; (b) entropy S; 
(c) specific heat C v under constant volume 



calculated free energy, entropy, specific heat under con- 
stant volume, and Debye temperature via the new classi- 
cal force field are in excellent agreement with those from 
ab initio calculations. 

Overall, our results indicate that the force field pre- 
sented provides a very accurate description of atomic in- 
teractions in rutile TiC>2. The force field can be applied 
with some confidence to studies of bulk properties of ti- 
tania. While further testing is necessary, we can also 
hope that it provides an improved description of surface 
properties. If this is the case, it can be used to study, 
not only surfaces, but also nanocrystals and nanowires, 
all of which have huge technological importance for their 
catalytic and optical properties. 
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We stress that no experimental data, other than the 
rutilc crystal symmetry, has been used to parameterize 
this potential, therefore, we can be confident that its pre- 
dictive capability is based on a robust description of the 
ions' potential-energy surface. 
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